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SECULAR EVOLUTION OF HD 12661: A SYSTEM CAUGHT AT AN UNLIKELY TIME 
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ABSTRACT 

The eccentricity evolution of multiple planet systems can provide valuable constraints on planet 
formation models. Unfortunately, the inevitable uncertainties in the current orbital elements can 
lead to significant ambiguities in the nature of the secular evolution. Integrating any single set of 
orbital elements inadequately describes the full range of secular evolutions consistent with current 
observations. Thus, we combine radial velocity observations of HD 12661 with Markov Chain Monte 
Carlo sampling to generate ensembles of initial conditions for direct n-body integrations. We find that 
any mean motion resonances are quite weak and do not significantly impact the secular evolution, 
and that current observations indicate circulation or large amplitude libration of the periapses. The 
eccentricity of the outer planet undergoes large oscillations for nearly all of the allowed two-planet 
orbital solutions. This type of secular evolution would arise if planet c had been impulsively perturbed, 
perhaps due to strong scattering of an additional planet that was subsequently accreted onto the star. 
Finally, we note that the secular evolution implied by the current orbital configuration implies that 
planet c spends ~ 96% of the time following an orbit more eccentric than that presently observed. 
Either this system is being observed during a relatively rare state, or additional planets are affecting 
the observed radial velocities and/or the system's secular eccentricity evolution. 
Subject headings: celestial mechanics — stars: individual (HD 12661) — planetary systems: formation 
— methods; n-body simulations, statistical 



1. INTRODUCTION 

The secular evolution of multi-planet systems has 
become a topic of considerable imp ortance for con- 
straining planetary forma ti on theories (iFord et al.ll2005 ; 
Adams fc LauEhlinI 120061 iBarnes fc Greenbergl l2006bl : 
Sandor et a l. 2007). We investigate the secular evo- 



lution of the two giant planets orbiting HD 12661, a 
~ 1.136M0 star. Planet HD 12661 b has a semima- 
jor axis of ~ 0.83AU and a moderate eccentricity, and 
planet HD 12661 c follow s a nearly circular orbit about 
three times further away ({Fischer et al. 2003). This sys- 
tem has inspired several dynamical studies. A brief sur- 
vey of their results demonstrates the importance of con- 
sider ing the uncertainty in th e curre nt orbital configura- 
tion. iKiseleva-Eggleton et al.l ()2002l ) considered coplanar 
edge-on systems and found that the then-current best-fit 
solution was chaotic, behavior that is likely due to close 
proxi mity to the 9:2 m ean- motion resonance (M MR) . 
Both lGozdziewskil (|2003f ) and lLee fc Peal^ (|2003f ) found 
that the system was close to the 11:2 MMR and that 
the periapses underwent large amplitude libration for 
the edge -on case, as w ell as for a broad range of incli- 
nations. I Zhou fc SunI (2003. ) claimed that aligned con- 
figurations were le ss likely to b e chao tic and thus more 
likely to be stable. iGozdziewskil ()2003[ ) found the system 
could be chaotic, but still stable for ~Gyr thanks to a 
secular apsida l lock about an anti aligned configuration. 
[Gozdziewski fc Macieiewskil (l2003fl performed an inde- 
pendent analysis of the lFischer et al.l (|2003f l radial veloci- 
ties to find a better orbital fit that placed the system very 
near (and likely in) the 6:1 MMR that could be stabilized 
by secular apsidal lock for a broad range of inclinations. 
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Despit e the different orbital solutions, bo t h iLee fc Peald 
(|2003f) and iGozdziewski fc Macieiewskil (l2003() found 
that the periapses were more likely to librate about an 
antialigned configuration, but that librations about an 
aligned configuration were also possible, particularly for 
systems with large inclinations. However. iJi et al.l (|2003l ) 
reported that aligned and antialigne d configuratio n s were 
nearly equally likely. Subsequently, iButler et al.l ()2006f ) 
published an orbital solut ion with a ratio of or b ital pe- 
riods approaching 13:2. iBarnes fc Greenbergl (|2006af l 
found that this orbital solution places the system very 
near the the boundary between librating and circulating 
modes of secular evolution. The revised orbital solutions 
for HD 12661 b fc c and the lack of consensus regarding 
the system's secular evolution both motivate an updated 
dynamical study. Further, the historical range of orbital 
configurations illustrates the importance of properly ac- 
c ounting for uncertain ties in orbital determinations. 

IGozdziewskil (2003) found that the classical secular 
theory gave only a crude approximation to the secular 
evolu tion due to large ec c entric ities a nd the 11:2 MMR. 
Both IVeras fc Armitagi (|2007t ) and iLibert fc HenrardI 
(j2007l ) caution against using low-o rder secular t heorie s 
to model the system' s behav ior. Lee fc Peale (2003'), 
iRodriguez fc Gallardol (|2005f ). and Libert fc HcnrarJ 
(I2007D ibund that either the octupole or high-order 
Laplace-Lagrange approximation described the secular 
evolution of HD 12661 well and that the secular dynam- 
ics is not significantly affected by the 5:1, 11:2, and 6:1 
MMRs. Although we do not expect near-resonant terms 
to significantly affect the secular evolution, we use direct 
n-body integrations to account for the full dynamics, in- 
cluding all possible MMRs. 

We determine the mode(s) of secular evolution consis- 
tent with observations and the requirement of long-term 
stability. Our study improves upon previous studies by 
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combining a Bayesian analysis of an updated set of ra- 
dial velocity observations with direct n-body integrations 
to examine HD 12661's stability and secular evolution. 
We characterize the full range of orbital histories that 
are consistent with observations and reject solutions that 
do not exhibit long-term orbital stability. Thus, we can 
make quantitative statements about the relative proba- 
bility of different modes of secular evolution. 

2. METHODS 

We reanalyze an updated set of radial velocity ob- 
servations (Wright et al. 2008) assuming that the stel- 
lar velocity is the superposition of two Keplerian or- 
bit plus noise. First, we perform a brute force search 
over parameter space to verify that there are no qualita- 
tively different and comparably likely solutions. Work- 
ing in a Bayesian framework, we generate an ensemble of 
~ 5 X 10^ orbital solutions usi ng Markov Cha in Monte 
Carlo (MCMC; Ford 2005, 2006; Gregorvl[2007allbl 'l . Each 
state of the Markov chain includes the orbital period 
(P), velocity amplitude (K), eccentricity (e), argument 
of pericenter measured from the plane of the sky (w) , and 
mean anomaly at a given epoch (u) for each planet. We 
use the priors, candidate transition functions, auto mated 
step s ize control, and convergence tests described in lFordI 
()2006t ). In particular, we assume a prior for the velocity 
"jitter" (aj) oi p{a.j) oc [I + aj/iluis'^)]-^. 

We randomly select subsamples of orbital solutions to 
be investigated with long-term n-body integrations and 
consider a range of possible inclinations (i) and ascend- 
ing nodes (Q). We generate each orbit from an isotropic 
distribution (i.e., uniform in cosi and f2) and divide the 
resulting systems into bins based on the relative incli- 
nation between the two orbital planes (irei)- The bins 
are 1) irei = 0° (coplanar), 2) 0° < i^d < 30°, 3) 30° < 
irel < 60°, 4) 60° < irel < 90°, 5) 90° < irel < 120°, 6) 
120° < irci < 150°, and 7) 150° < irei < 180°. We inte- 
grate 2,000 systems in each of the first four bins and 
500 in each of the last three bins. We can reweight 
our simulations so as to determine the probability of 
the different modes of secular evolution for various as- 
sumed distributions of relative inclinations. From each 
set {P,K,e,uj,i,fl,u), we generate the planet mass (m) 
and semi-maior axi s (a) using a Jacobi coordinate system 
(|Lee fc Pealei[2003h . 

We used the h ybrid symplectic integrator of Mercury 
(jChambersI |1999[ ) to integrate each set of initial condi- 
tions for at least 1 Myr. Based on a smaller series of 
10 Myr integrations, we found that the vast majority of 
instabilities are manifest within 1 Myr. We classified sys- 
tems as "unstable" if, for either planet, Cmax — a,min > 
TQi, where Cmax and Omin represent the maximum and 
minimum values of the semimajor axis, is the ini- 
tial value of the semimajor axis, and r — 0.3. We dis- 
carded each set of initial conditions that was found to be 
unstable, and analyzed the properties of the remaining 
"stable" systems. Although a small fraction of our "sta- 
ble" systems might exhibit instability if integrated for 
much longer timescales, our criteria avoids miscategoriz- 
ing a system exhibiting bounded chaos fe.g.. lGozdziewskil 
|2003() as unstable. We manually verified that the above 
criteria gives reasonable results and that various choices 
of T G [0.10,0.30] make just a few percent difference in 
the number of simulations labeled as stable. The percent 



of stable systems in our 7 bins are: 1) 100%, 2) 99.7%, 3) 
93.3%, 4) 16.4%, 5) 0.2%, 6) 16.7%, and 7) 98.6%. We 
independently affirmed the trend exhibited by these sta- 
bility percentages by performing a smaller, additional set 
of simulations by using an ensemble of initial conditions 
generated from a different MCMC code that accounted 
for planet-planet interactions with a Hermite integrator. 

3. RESULTS 



Several previous studi e s of HD 12661 (e .g. iJi et al.l 
[200l iLee fc Peald [200l iZhou fc SunI 120031) and other 



planetary systems Ichiang et all l2001t iMalhotral l2002t 
iFord et al.|[2005t iBarnes fc Greenberdl2006b[ ) have high- 
lighted the importance of the apsidal angle (Azu = 
rif, -|- tjf, — — LOc), which is useful for describing the 
secular dynamics of systems with small relative inclina- 
tions. Since we consider systems with a wide range of 
relative inclinations, we instead focus on Aw', the angle 
between the periastron directions projected onto the in- 
variable plane. Typically, the apsidal angle is classified as 
circulating, librating about 0°, or librating about 180°. 
By inspecting plots of apsidal angle evolution for indi- 
vidual systems, we find that some systems spend most 
of the time librating about one center, but occasionally 
the apsidal angle circulates for a short period of time. 
Therefore, we consider two summary statistics describ- 
ing the secular evolution: the root mean square (RMS) 
of Atn' and the mean absolute deviation (MAD) of Aro' 
about the libration center, where the "libration center" 
is defined as the angle about which the RMS Anj' is min- 
imized. For a system undergoing small amplitude libra- 
tion about either center, the MAD and RMS deviation 
are of order the libration amplitude. For a system under- 
going uniform circulation, the MAD approaches 90° and 
the RMS deviation approaches 103.92°. These statistics 
can be used to identify the mode of secular evolution in 
clear cut cases and provide a quantitative measure that 
is well-defined even for systems with complex evolution. 

We calculate both measures using orbital elements 
measured every 1000 yr. Although the libration am- 
plitude (maximum absolute deviation of Aw') can be 
sensitive to the frequency of sampling and leads to ambi- 
guities, both RMS Aro' and MAD Aw' are more robust 
statistics that allows us to describe the secular evolu- 
tion with a simple quantitative measure, regardless of 
whether the system is librating, circulating, or switching 
between regimes due to short-term perturbations. This 
robustness is particularly important for studying systems 
where there is little distinction between the librating and 
circulating regimes, due to one planet's orbit periodically 
becoming nearly circular. 

For each n-body integration, we determine the libra- 
tion center and calculate both the RMS and MAD of 
Aw' about that center. We discover that the RMS 
Aw' ranges from 68° — 85° and MAD Aw' ranges from 
63°— 82°, depending on the relative inclination (see Table 
[T]). The libration centers are preferentially anti-aligncd 
for prograde, nearly-coplanar systems, and transition to 
almost entirely aligned for near-coplanar retrograde sys- 
tems. The stability is highly dependent on initial relative 
inclination; very few highly inclined retrograde systems 
were stable, indicating that non-secular perturbations in- 
fluence the long-term dynamics for some relative inclina- 
tions. Although large values of RMS Aw' and MAD Aw' 
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Fig. 1. — Values of ci (dots), C2 (triangles), C3 (squares), and 
C4 (crosses) and their standard deviations, for coplanar systems, 
binned according to the sine of the inclination of the invariable 
plane. 



may be consistent with non-uniform circulation, visual 
spot checks of individual systems all suggest libration. 

If a system of two planets on nearly circular orbits is 
excited impulsively, then the secular evolution will fol- 
low the boundary between the circulating and librating 
regimes and cause o ne orbit to repe atedly return to a 
nearly circular orbit (jMalhotrall2002f ). In some systems, 
this secular evol ution can be used to constrain the sys- 
tem's fo rmation (iFord et al.ll2005D. iBarnes fc Greenberd 
(|2006bf) found that the iButler et all (|2006[) orbital solu- 
tion implies the eccentricity of HD 12661 c comes par- 
ticularly close to zero, as pa rameterized by their e pa- 
rameter (Eq. 1 of iBarnes fc Grecnberg 2006b). Because 
they considered only the published orbital solution in an 
edge-on, coplanar orientation, they were unable to asses 
the finding's robustness to uncertainties in the orbital pa- 
rameters or non-edge-on systems. In order to determine 
what fraction of stable systems consistent with observa- 
tions result in one planet returning to a nearly circular 
orbit, we calculate four statistics from each of our n-body 
integrations. The first two are ci = emin.c/ emax,c and 
C2 = emin,b/smax.bi which represent the ratio of the mini- 
mum to maximum eccentricity for each planet. The third 

(C3) is equal t o [2min(efeee)] /(Xmax-Xrmn+ ymax-ymin). 

as defined by IBarnes &: Greenbergl ()2006bl ) , where x = 
CfjCc sin (Aro) and y = e^ec cos (Aw). The fourth is 
C4 = min(ef,ec)/max(ef,ec), which is a similar measure 
that includes the eccentricities of both planets in a ro- 
tationally symmetric manner that is independent of the 
systems' orientation (unlike C3). 

Figure [T] plots the mean and standard deviation of 
each measure (ci - dots, C2 - triangles, C3 - squares, C4 - 
crosses) for coplanar systems which have been binned 
according to the sine of the inclination of the invari- 
able plane. The values ci, C3 and C4 all maintain values 
< 0.1 for all bins of relative inclination for prograde sys- 
tems, and < 0.15 for retrograde systems. The large val- 
ues of C2 (triangles) imply that HD 12661 b maintains a 
moderately eccentric orbit throughout the planet's evo- 
lution. Regardless, the measure C4 demonstrates that 
the system lies near the boundary of libration and cir- 



FlG. 2. — Fraction of time for which each planet's eccentricity 
exceeded a given threshold (x-axis) averaged over allowed orbital 
solutions. The curves with dots show results for the coplanar case, 
but assuming an isotropic distribution of inclinations relative to 
the plane of the sky. The dashed vertical lines indicate the me- 
dian eccentricity values (0.031 and 0.361) for the planets' initial 
eccentricities, and their 5%-95% uncertainty ranges are indicated 
by the horizontal solid segments. HD 12661 c's orbital eccentricity 
exceeds its current value for approximately 96% of its evolution. 

culation. For the vast majority of systems, C3 is over 
one order of magnitude gre ater than 0.003 (as found by 
IBarnes fc Greenberel[2006b[) . 

The finding that planet c undergoes large eccentric- 
ity oscillations is remarkable since it currently has an 
eccentricity of only ~ 0.02. We investigate this observa- 
tion further by plotting the fraction of simulation time 
that each planet's eccentricity exceeded a given value for 
the coplanar systems (Fig. [2]). The dotted vertical lines 
indicate representative values for the planets' initial ec- 
centricities. HD 12661 c's orbital eccentricity exceeds 
its current value for approximately 96% of its evolution. 
Thus, the near-circular orbit we now observe is a rare 
occurrence. 

4. CONCLUSION 

Several previous studies had suggested HD 12661 b & c 
are in a secular apsidal lock and are undergoing large am- 
plitude librations. However, previously, even the orbital 
period of HD 12661 c was poorly determined. Now that 
radial velocity observations span over two orbital peri- 
ods of HD 12661 c, the current orbital period of planet 
c is well-constrained. An exhaustive search for libration 
of all possible 11:2, 6:1, and 13:2 resonant angles found 
that in only a few percent of all simulations performed, 
large amplitude libration lasts over several 10^ yr. The 
current orbital elements of the two giant planets imply 
that the system lies near the boundary of circulating and 
librating regimes, regardless of the unknown orbital incli- 
nations. This behavior causes the eccentricity of planet 
c to spend most of its time with a significant eccentricity 
(cc > 0.2). By computing the projected apsidal angle on 
the invariable plane, we find that the RMS Aw' ranges 
from 68° - 85° and MAD Aw' ranges from 63° - 82°, 
for any stable initial relative inclination, and planet c's 
eccentricity would be greater than its current value 96% 
of the time. 

The most straightforward interpretation is that the 
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eccentricity of HD 12661 b was excited by an impul- 
sive perturbation and planet c's eccentricity is peri- 
odically excited via secular perturbations, similar to 
the history of the v And c & d system ([Malhotral 
I2OO2I: iFord et all l2005| ). Such an impulse could be 
the result of planet-planet scat tering involving at 
least one additional planet (e .g. iRasio fc FordI 119961 : 
IWeidenschilling fc Marza7illl996D . The current eccentric- 
ity of planet b suggests that the putative scattered planet 
would have had a mass roughly half that of planet b 
(jFord fc Rasiol[2008^ . If the additional planet had been 
scattered outwards, then there is a significant chance that 
it would have altered the secular evolution of planet c 
(jBarnes fc Gree nbcrg 2006b). Given the semi-major axis 
of planet b, it is possible that the scattered planet was 
accreted onto the star. If the instability occurred after 
the star had reached the main sequence, then this in- 
stability could contribute to the star's large surface high 
metalicity. One potential concern for this mechanism is 
whether the timescale for the apoastron of the scattered 
planet to be lowered would be significantly shorter than 
the ~ 10** — lO^yr timescale for secular evolution of plan- 
ets b & c. 

There is an alternative interpretation of our results. 
One might think that we are unlikely to observe a sys- 
tem at such a rare time, casting doubt upon the orbital 
fit and thus the inferred secular evolution. One would 
expect that one of the 30 currently observed multiple 
planet systems is transiently passing through such an un- 
usual state. However, there are several ways in which a 
system could be "unusual" so that the standard caveats 
of a posteriori statistics apply. Given the current ra- 
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dial velocity observations, there are few alternatives. We 
consider it unlikely that our global search missed a qual- 
itatively different pair of orbits due to our assumption of 
non-interacting Keplerian orbits. More plausibly, one or 
more undetected planets could be affecting the orbital so- 
lution and/or the secular evolution of the system. After 
accounting for planets b fc c, we find no statistically sig- 
nificant periodicities in the current radial velocity obser- 
vations. Current observations are consistent with a small 
long-term acceleration, but zero acceleration falls within 
the 95% credible interval. Further, models with a long- 
term acceleration do not result in significantly different 
orbital parameters. Our Bayesian analysis assuming two 
planets and uncorrelated Gaussian noise constrains the 
jitter to be ~ 3.4 ± 0.7m/s. Further radial velocity ob- 
servations can test whether a portion of this jitter is due 
to additional planets. 
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TABLE 1 
Results of N-body Integrations 





(RMS) aligned fraction 


(MAD) aligned fraction 


RMS Aro' 


MAD Aro' 


Stable 


Isotropic 


0° 


35.9% 


31.7% 


84.9° ±8.5° 


81.1° ± 7.8° 


100% 




0° - 30° 


39.3% 


35.2% 


83.8° ±9.4° 


79.5° ± 8.6° 


99.7% 


6.7% 


30° - 60° 


56.8% 


51.8% 


84.2° ± 10.4° 


79.7° ±9.8° 


93.3% 


18.3% 


60° - 90° 


56.0% 


52.0% 


83.8° ± 10.5° 


79.6° ±9.9° 


16.4% 


25.0% 


90° - 120° 










0.2% 


25.0% 


120° - 150° 


95.2% 


95.2% 


74.7° ± 12.1° 


69.3° ± 11.1° 


16.7% 


18.3% 


150° - 180° 


99.6% 


99.2% 


68.4° ± 12.8° 


63.0° ± 12.1° 


98.6% 


6.7% 


Prograde, isotropic 
Retrograde, isotropic 










69.8% 
38.5% 


80.0% 
50.0% 



Note. — Columns 2 and 3 each list the RMS and MAD percent of stable systems with a libration center of 0° instead 
of 180° for the range of relative inclinations in column 1. Columns 4 and 5 report the RMS and MAD variations 
for systems which librate around the most common center for that range of initial relative inclinations. Column 6 
lists the percent of stable systems, and Column 7 hsts the percent of systems with the given range of inclination for 
an isotropic distribution. The bottom two rows present aggregated stability statistics for the isotropic prograde and 
retrograde distributions of orbits. 2000 simulations were performed for each prograde relative inclination bin, and 
500 simulations for each retrograde bin. 



